NetCoMi

For more details on using the NetCoMi package see https://github.com/stefpeschel/NetCoMi.

## Warning: package 'phyloseq' was built under R version 3.6.1
## [1] '1.30.0'

Networks comparison

The following code was run on a cluster.

# # load phyloseq data
# load('/n/home05/ajsommer/NetCoMi_cluster/ps_to_net.RData')

# gut_split <- metagMisc::phyloseq_sep_variable(ps_prune, "W")

# net_W <- netConstruct(gut_split$`0`, gut_split$`1`, verbose = 2,
# filtTax = "highestVar",
# filtTaxPar = list(highestVar = 50),
#                            measure = "spieceasi",
#                           measurePar = list(method = "glasso",
#                                             nlambda=20,
#                                             pulsar.params=list(rep.num=50)),
#                           normMethod = "none", zeroMethod = "none",
#                           sparsMethod = "none", seed = 123456, matchDesign = c(1,1))


# props_W <- netAnalyze(net_W, clustMethod = "cluster_fast_greedy")

We plot the results of the netAnalyze function.

phyl_ps_prune <- as.factor(tax_table(ps_prune)[, 'Phylum'])
names(phyl_ps_prune) <- rownames(tax_table(ps_prune)[, 'Phylum'])
phylcol <- rainbow(length(unique(phyl_ps_prune)))

set.seed(1)
props_W_igraph <- graph_from_adjacency_matrix(props_W$input$assoMat1, 
                                              #mode = "undirected",
                                              weighted = TRUE)
FA_coord <- layout.forceatlas2(props_W_igraph, iterations=2000, plotstep=2001)
rownames(FA_coord) <- names(phyl_ps_prune)

plot(props_W, 
     layout = FA_coord,
     # layout = "layout_with_fr",
     sameLayout = TRUE, layoutGroup = 1, labels = FALSE, 
     featVecCol = phyl_ps_prune, nodeColor = "feature",nodeTransp = 55,
     borderCol = "gray40", highlightHubs = FALSE,nodeSize = "normCounts", 
     nodeSizeSpread = 1, cexNodes = 5, cexTitle = .8,
     edgeTranspLow = 70, edgeTranspHigh = 50,
     mar = c(3,5,3,2), groupNames = c("Polluted air", "Clean air"))
legend("center", inset = c(-1,-10), 
       legend=levels(phyl_ps_prune),
       col=phylcol, bty="n", text.col = phylcol, cex = .5) 
legend("bottom", inset = c(0,-2),
       legend=c("positive","negative"), lty = 1,
       col=c("darkgreen","red"), bty="n", cex =.6)

Modularity

Vizualisation

nclust <- max(as.numeric(names(table(props_W$clustering$clust1))))+1
col <- polychrome(nclust)
names(col) <- 1:nclust

clusters <- as.factor(c(props_W$clustering$clust1,props_W$clustering$clust2))

col_assign <- col[clusters]

names(col_assign) <- names(clusters)

# which(tax_table(ps_prune)[,"Genus"] == "Blautia")
# tax_table(ps_prune)[which(tax_table(ps_prune)[,"Genus"] == "Bacteroides"),]
# 248

shapeVec <- rep(1, ncol(otu_table(ps_prune)))
shapeVec[248] <- 2
names(shapeVec) <- names(phyl_ps_prune)

plot(props_W,
     layout = FA_coord,
     # layout = "layout_with_fr",
     sameLayout = TRUE, layoutGroup = 1, 
     labels = list(props_W$clustering$clust1, props_W$clustering$clust2), 
     labelFont = 1, cexLabels = 2,
     nodeColor = "cluster",
     colorVec = col,
     nodeTransp = 40, 
     borderCol = "gray40", highlightHubs = FALSE,
     nodeSize = "fix", cexNodes = 3, 
     nodeShape = c("circle", "triangle"),
     featVecShape = shapeVec,
     edgeTranspLow = 80, edgeTranspHigh = 50,
     groupNames = c("Polluted air", "Clean air"))

Shared nodes among modules

dat_graph <- data.frame(id = names(props_W$clustering$clust1), 
                        com_p = props_W$clustering$clust1,
                        com_c = props_W$clustering$clust2)

dat_graph[,2] <- paste0(dat_graph[,2], "_P")
dat_graph[,3] <- paste0(dat_graph[,3], "_C")

inc_com_sum <- ddply(data.frame(dat_graph), .(com_p, com_c), 
                     summarise, sum = length(id))

g_sum <- graph.data.frame(inc_com_sum[,c(1,2)], directed = F)
V(g_sum)$type <- V(g_sum)$name %in% inc_com_sum[,1] #the first column of edges is TRUE type

# node size represents the number of OTUs in each module
table(props_W$clustering$clust1)
## 
##  0  1  2  3  4  5  6  7  8  9 10 11 12 
## 26 62 94 86 55 25 45 98  4 12  3  2  3
table(props_W$clustering$clust2)
## 
##   0   1   2   3   4   5   6   7   8 
##  10  65  68  76 113 104  30  47   2
# add community name (without T or C)
V(g_sum)$com_nr <- as.numeric(unlist(strsplit(V(g_sum)$name, "_"))[c(TRUE,FALSE)])
V(g_sum)$com_size[V(g_sum)$type == TRUE] <- table(props_W$clustering$clust1)[V(g_sum)$com_nr[V(g_sum)$type == TRUE]] # control sizes
## Warning in V(g_sum)$com_size[V(g_sum)$type == TRUE] <-
## table(props_W$clustering$clust1)[V(g_sum)$com_nr[V(g_sum)$type == : number of
## items to replace is not a multiple of replacement length
V(g_sum)$com_size[V(g_sum)$type == FALSE] <- table(props_W$clustering$clust2)[V(g_sum)$com_nr[V(g_sum)$type == FALSE]] # treated sizes
## Warning in V(g_sum)$com_size[V(g_sum)$type == FALSE] <-
## table(props_W$clustering$clust2)[V(g_sum)$com_nr[V(g_sum)$type == : number of
## items to replace is not a multiple of replacement length
V(g_sum)$color <- V(g_sum)$type
V(g_sum)$color=gsub("FALSE","#1CFFCE",V(g_sum)$color)
V(g_sum)$color=gsub("TRUE","#66B0FF",V(g_sum)$color)
E(g_sum)$weight <- as.numeric(inc_com_sum[,3])

layout_bi_rows <- layout_as_bipartite(g_sum, hgap = 50, vgap = 60, maxiter = 1000)

legend_cats <- data.frame(size_circle = V(g_sum)$com_size/3,
                          size_text = V(g_sum)$com_size)
legend_cats <- legend_cats[order(legend_cats$size_circle), c(1,2)]
# dim(legend_cats)
legend_cats <- legend_cats[c(1,5,10,15,22),]

# order by vertex size
par(mfrow=c(1,1))
plot(g_sum, vertex.size = V(g_sum)$com_size/3, edge.width=E(g_sum)$weight/5, 
     layout=layout_bi_rows[,2:1], vertex.label=V(g_sum)$com_nr,
     vertex.label.dist=5, vertex.label.degree = pi*V(g_sum)$type, 
     margin = -.2, asp = 2, edge.arrow.size=.7)

legend(x=-1.5, y=-1.1, c("W = 0 (High)", "W = 1 (Low)"), pch=21,
       col="#777777", pt.bg=c("#66B0FF","#1CFFCE"), pt.cex=2, cex=.8, bty="n", ncol=1)

legend(x=-.03, y=-.9, legend=legend_cats$size_text, pt.cex=legend_cats$size_circle/9, bty="n", cex = .8,
       col="#777777", pch=21, pt.bg="lightgray", ncol = 2, y.intersp = 1.7, x.intersp = 1.5, 
       title = "Module size (nb. of nodes):")

Module 7 split

names_clust <- as.character(names(which(props_W$clustering$clust1 == 7)))
labels_phyl <- substr(as.character(phyl_ps_prune[names_clust]), 1, 3)

plot(props_W, 
     # layout = FA_coord,
     sameLayout = TRUE, layoutGroup = 1, 
     labels = labels_phyl, 
     labelFont = 1, cexLabels = 1.5,
     featVecCol = phyl_ps_prune[names_clust], 
     nodeColor = "cluster",
     colorVec = col,
     nodeTransp = 40, 
     borderCol = "gray40", highlightHubs = FALSE,
     nodeSize = "fix", cexNodes = 1, 
     nodeFilter = "names",
     nodeShape = c("circle", "triangle"),
     featVecShape = shapeVec,
     nodeFilterPar = names_clust,
     edgeTranspLow = 80, edgeTranspHigh = 50,
     groupNames = c("Polluted air", "Clean air"))

Correlation structure of module 7

col1 <- colorRampPalette(c("#7F0000", "red", "#FF7F00", "yellow", "white",
                           "cyan", "#007FFF", "blue", "#00007F"))

cor_0 <- props_W$input$assoMat1[names_clust,names_clust]
cor_1 <- props_W$input$assoMat2[names_clust,names_clust]

par(mfrow = c(1,2))
corrplot(cor_0, method = "color", tl.col = "grey", mar=c(0,0,5,0),
         col = col1(50), title = "Polluted air")
corrplot(cor_1, method = "color", tl.col = "grey", mar=c(0,0,5,0),
         col = col1(100), title = "Clean air")